# replication file for
# "Ideological Positions and Committee Chair Appointments"
# published in Legislative Studies Quarterly
# authors: Jochen Rehmert & Naofumi Fujimura
 
# replicates
# Figure 1

# needs:
# ChaiR_stata2.dta
# coef_dist_pc2_re_dist_ldp.csv
# coef_hp_pc2_re_dist_ldp.csv
# coef_oth_pc2_re_dist_ldp.csv
# vcov_dist_pc2_re_dist_ldp.csv
# vcov_hp_pc2_re_dist_ldp.csv
# vcov_oth_pc2_re_dist_ldp.csv


# directory 
wd <- "C:/..."
setwd(wd)
# prepare data
library(readstata13)
dat <- readstata13::read.dta13("ChaiR_stata2.dta")
dat <- dat[dat$imp == 0,]

# generate covariates 
dat$ln_pop <- log(dat$pop_density)
dat$zombie <- ifelse(dat$result == 2, 1,0)
dat$exp_high <- ifelse(dat$cab_exp >= 1 | dat$jud_exp >= 1 | dat$security_exp >= 1 | dat$fin_exp  >= 1 | dat$foreign_exp >= 1 | dat$basic_exp >= 1 | dat$health_exp >= 1
                       , 1,0)
dat$exp_foreign <- ifelse(dat$cab_exp >= 1 | dat$jud_exp >= 1 | dat$security_exp >= 1 |  dat$foreign_exp >= 1 , 1,0)

# list files in folder
files.coef <- list.files(pattern = "coef")
files.vcov <- list.files(pattern = "vcov")

# function to read coefficients & variance matrices from stata
read.files <- function(vars = c("all","hp","fa","dist", "re_dist_ldp")){
  
  
  tmp.1 <- read.csv(files.coef[which(grepl(vars, files.coef))][3], sep = ",") # other chair
  tmp.2 <- read.csv(files.coef[which(grepl(vars, files.coef))][2], sep = ",") # high chair
  tmp.3 <- read.csv(files.coef[which(grepl(vars, files.coef))][1], sep = ",") # dist chairs
  
  coefs.names1 <- as.character(tmp.1[seq(from = 3, to = nrow(tmp.1), by = 2), 1])
  coefs.names1 <- as.character(gsub("=","",coefs.names1))[-length(coefs.names1)]
  coefs.names2 <- as.character(tmp.2[seq(from = 3, to = nrow(tmp.2), by = 2), 1])
  coefs.names2 <- as.character(gsub("=","",coefs.names2))[-length(coefs.names2)]
  coefs.names3 <- as.character(tmp.3[seq(from = 3, to = nrow(tmp.3), by = 2), 1])
  coefs.names3 <- as.character(gsub("=","",coefs.names3))[-length(coefs.names3)]
   
  coefs.est1 <- as.character(tmp.1[seq(from = 3, to = nrow(tmp.1), by = 2), 2])
  coefs.est1 <- as.numeric(as.character(gsub("=","",coefs.est1)))[-length(coefs.est1)]
  coefs.est2 <- as.character(tmp.2[seq(from = 3, to = nrow(tmp.2), by = 2), 2])
  coefs.est2 <- as.numeric(as.character(gsub("=","",coefs.est2)))[-length(coefs.est2)]
  coefs.est3 <- as.character(tmp.3[seq(from = 3, to = nrow(tmp.3), by = 2), 2])
  coefs.est3 <- as.numeric(as.character(gsub("=","",coefs.est3)))[-length(coefs.est3)]
   
  names(coefs.est1) <- coefs.names1
  names(coefs.est2) <- coefs.names2
  names(coefs.est3) <- coefs.names3
  
  
  tmp.vcov1 <- read.csv(files.vcov[which(grepl(vars, files.vcov))][3], sep = ",")[-1,] # other chair
  tmp.vcov2 <- read.csv(files.vcov[which(grepl(vars, files.vcov))][2], sep = ",")[-1,] # high chair
  tmp.vcov3 <- read.csv(files.vcov[which(grepl(vars, files.vcov))][1], sep = ",")[-1,] # dist chairs
  
  
  tmp.vcov1.tmp <- tmp.vcov1[-c(1:2),-c(1:2)]
  tmp.vcov1.tmp <- tmp.vcov1.tmp[-which(row.names(tmp.vcov1.tmp) == "24"),]
  tmp.vcov2.tmp <- tmp.vcov2[-c(1:2),-c(1:2)]
  tmp.vcov2.tmp <- tmp.vcov2.tmp[-which(row.names(tmp.vcov2.tmp) == "24"),]
  tmp.vcov3.tmp <- tmp.vcov3[-c(1:2),-c(1:2)]
  tmp.vcov3.tmp <- tmp.vcov3.tmp[-which(row.names(tmp.vcov3.tmp) == "24"),]
 
  var.names1 <- as.character(tmp.vcov1[,1])[-c(1:2)]
  var.names1 <- var.names1[-21]
  
  var.names2 <- as.character(tmp.vcov2[,1])[-c(1:2)]
  var.names2 <- var.names2[-21]
  var.names3 <- as.character(tmp.vcov3[,1])[-c(1:2)]
  var.names3 <- var.names3[-21]
  
  colnames(tmp.vcov1.tmp) <- row.names(tmp.vcov1.tmp) <- var.names1
  colnames(tmp.vcov2.tmp) <- row.names(tmp.vcov2.tmp) <- var.names2
  colnames(tmp.vcov3.tmp) <- row.names(tmp.vcov3.tmp) <- var.names3
  
  tmp.vcov1 <- data.frame(tmp.vcov1.tmp)
  tmp.vcov2 <- data.frame(tmp.vcov2.tmp)
  tmp.vcov3 <- data.frame(tmp.vcov3.tmp)
  
  for(cc in 1:ncol(tmp.vcov1)){tmp.vcov1[,cc] <- as.numeric(as.character(tmp.vcov1[,cc]))}
  for(cc in 1:ncol(tmp.vcov2)){tmp.vcov2[,cc] <- as.numeric(as.character(tmp.vcov2[,cc]))}
  for(cc in 1:ncol(tmp.vcov3)){tmp.vcov3[,cc] <- as.numeric(as.character(tmp.vcov3[,cc]))}
  
  return(list(coef.pc2.oth = coefs.est1,
              coef.pc2.high = coefs.est2, 
              coef.pc2.dist = coefs.est3,
              vcov.pc2.oth = tmp.vcov1,
              vcov.pc2.high = tmp.vcov2,
              vcov.pc2.dist = tmp.vcov3))
}



com.hp <- read.files("re_dist_ldp")

### predict probabilities ###
library(MASS)
seeds = 1234
# -- PC2: high policy -- #
set.seed(seeds)
n.sim = 1000
b.sim <- mvrnorm(n.sim, as.vector(com.hp$coef.pc2.high)[-c(1,((length(com.hp$coef.pc2.high)-1):length(com.hp$coef.pc2.high)))],
                 as.matrix(com.hp$vcov.pc2.high)[-nrow(com.hp$vcov.pc2.high), -ncol(com.hp$vcov.pc2.high)])
# create scenario with low value on ideological distance
scen.0 <- data.frame(cbind(0,0,0, 0, 0, 0, 0, 0, 0, 0,                      
                           female = median(dat$female, na.rm = TRUE),
                           age = mean(dat$age, na.rm = TRUE),
                           term = 5,
                           exp_high = 1,
                           zombie = median(dat$zombie, na.rm = TRUE),
                           pr_only = median(dat$pr_only, na.rm = TRUE),
                           ln_pop = mean(dat$ln_pop, na.rm = TRUE),
                           main_faction = 1,
                           dist_cab_pc2_ldp = 0,
                           1))
# create scenario with high value on ideological distance
scen.1 <- data.frame(cbind(0,0,0, 0, 0, 0, 0, 0, 0, 0,                      
                           female = median(dat$female, na.rm = TRUE),
                           age = mean(dat$age, na.rm = TRUE),
                           term = 5,
                           exp_high =1,
                           zombie = median(dat$zombie, na.rm = TRUE),
                           pr_only = median(dat$pr_only, na.rm = TRUE),
                           ln_pop = mean(dat$ln_pop, na.rm = TRUE),
                           main_faction = 1,
                           dist_cab_pc2_ldp = 1,
                           1))

# apply logit link function
lin.preds.0 <- as.matrix(scen.0) %*% t(b.sim)
preds.0 <- 1/(1+exp(-lin.preds.0))

lin.preds.1 <- as.matrix(scen.1) %*% t(b.sim)
preds.1 <- 1/(1+exp(-lin.preds.1))

# obtain average marginal effects
(cis.95.high <- apply(preds.1 - preds.0, 1, quantile, probs = c(0.025, 0.975))*100)
(cis.90.high <- apply(preds.1 - preds.0, 1, quantile, probs = c(0.05, 0.95))*100)
(ame.high <- apply(preds.1 - preds.0, 1, mean)*100)


# -- PC2: distributive committee -- *
set.seed(seeds)
b.sim.dist <- mvrnorm(n.sim, as.vector(com.hp$coef.pc2.dist)[-c(1,((length(com.hp$coef.pc2.dist)-1):length(com.hp$coef.pc2.dist)))],
                 as.matrix(com.hp$vcov.pc2.dist)[-nrow(com.hp$vcov.pc2.dist), -ncol(com.hp$vcov.pc2.dist)])

# create scenario with low value on ideological distance
scen.0 <- data.frame(cbind(0,0,0, 0, 0, 0, 0, 0, 0, 0,                      
                           female = median(dat$female, na.rm = TRUE),
                           age = mean(dat$age, na.rm = TRUE),
                           term = 5,
                           exp_high = 1,
                           zombie = median(dat$zombie, na.rm = TRUE),
                           pr_only = median(dat$pr_only, na.rm = TRUE),
                           ln_pop = mean(dat$ln_pop, na.rm = TRUE),
                           main_faction = 1,
                           dist_cab_pc2_ldp = 0,
                           1))
# create scenario with high value on ideological distance
scen.1 <- data.frame(cbind(0,0,0, 0, 0, 0, 0, 0, 0, 0,                      
                           female = median(dat$female, na.rm = TRUE),
                           age = mean(dat$age, na.rm = TRUE),
                           term = 5,
                           exp_high =1,
                           zombie = median(dat$zombie, na.rm = TRUE),
                           pr_only = median(dat$pr_only, na.rm = TRUE),
                           ln_pop = mean(dat$ln_pop, na.rm = TRUE),
                           main_faction = 1,
                           dist_cab_pc2_ldp = 1,
                           1))

# apply logit link function
lin.preds.0 <- as.matrix(scen.0) %*% t(b.sim.dist)
preds.0 <- 1/(1+exp(-lin.preds.0))

lin.preds.1 <- as.matrix(scen.1) %*% t(b.sim.dist)
preds.1 <- 1/(1+exp(-lin.preds.1))

# obtain average marginal effects
(cis.95.dist <- apply(preds.1 - preds.0, 1, quantile, probs = c(0.025, 0.975))*100)
(cis.90.dist <- apply(preds.1 - preds.0, 1, quantile, probs = c(0.05, 0.95))*100)
(ame.dist <- apply(preds.1 - preds.0, 1, mean)*100)

# -- PC2: other committee -- #
set.seed(seeds)
b.sim.oth <- mvrnorm(n.sim, as.vector(com.hp$coef.pc2.oth)[-c(1,((length(com.hp$coef.pc2.oth)-1):length(com.hp$coef.pc2.oth)))],
                 as.matrix(com.hp$vcov.pc2.oth)[-nrow(com.hp$vcov.pc2.oth), -ncol(com.hp$vcov.pc2.oth)])

# create scenario with low value on ideological distance
scen.0 <- data.frame(cbind(0,0,0, 0, 0, 0, 0, 0, 0, 0,                      
                           female = median(dat$female, na.rm = TRUE),
                           age = mean(dat$age, na.rm = TRUE),
                           term = 5,
                           exp_high = 1,
                           zombie = median(dat$zombie, na.rm = TRUE),
                           pr_only = median(dat$pr_only, na.rm = TRUE),
                           ln_pop = mean(dat$ln_pop, na.rm = TRUE),
                           main_faction = 1,
                           dist_cab_pc2_ldp = 0,
                           1))
# create scenario with high value on ideological distance
scen.1 <- data.frame(cbind(0,0,0, 0, 0, 0, 0, 0, 0, 0,                      
                           female = median(dat$female, na.rm = TRUE),
                           age = mean(dat$age, na.rm = TRUE),
                           term = 5,
                           exp_high =1,
                           zombie = median(dat$zombie, na.rm = TRUE),
                           pr_only = median(dat$pr_only, na.rm = TRUE),
                           ln_pop = mean(dat$ln_pop, na.rm = TRUE),
                           main_faction = 1,
                           dist_cab_pc2_ldp = 1,
                           1))

# apply logit link function
lin.preds.0 <- as.matrix(scen.0) %*% t(b.sim.oth)
preds.0 <- 1/(1+exp(-lin.preds.0))

lin.preds.1 <- as.matrix(scen.1) %*% t(b.sim.oth)
preds.1 <- 1/(1+exp(-lin.preds.1))

# obtain average marginal effects
(cis.95.oth <- apply(preds.1 - preds.0, 1, quantile, probs = c(0.025, 0.975))*100)
(cis.90.oth <- apply(preds.1 - preds.0, 1, quantile, probs = c(0.05, 0.95))*100)
(ame.oth <- apply(preds.1 - preds.0, 1, mean)*100)

# combine predicted probabilities in data.frame()
plotDat <- data.frame(rbind(cbind(ame = ame.high, committee ="High Policy Committee", low = cis.95.high[2] , high = cis.95.high[1], ci = 95),
                            cbind(ame = ame.high, committee ="High Policy Committee", low = cis.90.high[2] , high = cis.90.high[1], ci = 90),
                            cbind(ame = ame.dist, committee ="Distributive Committee", low = cis.95.dist[2] , high = cis.95.dist[1], ci = 95),
                            cbind(ame = ame.dist, committee ="Distributive Committee", low = cis.90.dist[2] , high = cis.90.dist[1], ci = 90),
                            cbind(ame = ame.oth, committee ="Other Committee", low = cis.95.oth[2] , high = cis.95.oth[1], ci = 95),
                            cbind(ame = ame.oth, committee ="Other Committee", low = cis.90.oth[2] , high = cis.90.oth[1], ci = 90)))

# make sure numeric variables are stored as numerics
plotDat$ame <- as.numeric(as.character(plotDat$ame))
plotDat$low <- as.numeric(as.character(plotDat$low))
plotDat$high <- as.numeric(as.character(plotDat$high))


##################################################################
# Figure 1: Ideological Distance from Executive Party Leadership #
##################################################################


par(mar = c(5,7,2,2) + 0.1)
plot(plotDat$ame[plotDat$ci == 95] , c(.9,1.1, 1.3), xlim = c(-7.2,2), ylim = c(0.8, 1.6) , pch = 19, xaxt = "n",
     xlab = "Change in Predicted Probabilities", yaxt ="n", frame = F, ylab = "")
segments(x0 = plotDat$low[plotDat$ci == 95], y0 = c(.9,1.1, 1.3), x1 =  plotDat$high[plotDat$ci == 95], y1 = c(.9,1.1, 1.3), lty = 2)
segments(x0 = plotDat$low[plotDat$ci == 90], y0 = c(.9,1.1, 1.3), x1 =  plotDat$high[plotDat$ci == 90], y1 = c(.9,1.1, 1.3))
lines(rep(0,2), c(0.7,1.4), col = "red", lty = 2)
axis(1, at = seq(-10,4,2), labels = seq(-10,4,2), col.ticks = "black", col = "white")
axis(2, at = c(1.5,.9,1.1, 1.3), labels = c("Type of  \nCommittee", "High Policy", "Distributive","Other"), col.ticks = "white", col = "white",
     las = 2)







